分位点回帰

自作関数での分位点回帰

# 損失関数1
loss.f1 <- function(tau = 0.5, y, a) {
    sum(ifelse((y - a) < 0, tau - 1, tau) * (y - a))
}

# 損失関数1の最小化で分位点を求める
y0 <- 0:10
opt <- optim(par = 0, fn = loss.f1, tau = 0.5, y = y0, method = "BFGS")
opt$par

## [1] 5

# 損失関数2：線形回帰
loss.f2 <- function(tau = 0.5, y, x, a) {
    sum(ifelse((y - a[1] - a[2] * x) < 0, tau - 1, tau) * (y - a[1] - a[2] * 
        x))
}

# 人工データ
set.seed(200)
x <- runif(100)
e <- rnorm(100, 0, 0.5)
y <- 5 + 2 * x + e
plot(x, y)

# 損失関数2をの最小化で分位点回帰
tau <- 0.5
opt <- optim(par = c(0, 0), fn = loss.f2, tau = tau, y = y, x = x, 
    method = "BFGS")
opt$par

## [1] 5.022 2.028

abline(a = opt$par[1], b = opt$par[2], col = "black")

tau <- 0.9
opt <- optim(par = c(0, 0), fn = loss.f2, tau = tau, y = y, x = x, 
    method = "BFGS")
opt$par

## [1] 5.683 1.965

abline(a = opt$par[1], b = opt$par[2], col = "red")

tau <- 0.1
opt <- optim(par = c(0, 0), fn = loss.f2, tau = tau, y = y, x = x, 
    method = "BFGS")
opt$par

## [1] 4.285 2.220

abline(a = opt$par[1], b = opt$par[2], col = "blue")

quantregパッケージによる分位点回帰

library(quantreg)

## Loading required package: SparseM

## Package SparseM (0.96) loaded.
##     To cite, see citation("SparseM")
## 

## 
## Attaching package: 'SparseM'
## 

## The following object(s) are masked from 'package:base':
## 
##     backsolve
## 

# 人工データ
set.seed(200)
x <- runif(100)
e <- rnorm(100, 0, 0.5)
y <- 5 + 2 * x + e
plot(x, y)

(fit <- rq(y ~ x, tau = 0.5))

## Call:
## rq(formula = y ~ x, tau = 0.5)
## 
## Coefficients:
## (Intercept)           x 
##       5.020       2.029 
## 
## Degrees of freedom: 100 total; 98 residual

abline(a = fit$coef[1], b = fit$coef[2], col = "black", lty = "dashed", 
    lwd = 3)

(fit <- rq(y ~ x, tau = 0.9))

## Call:
## rq(formula = y ~ x, tau = 0.9)
## 
## Coefficients:
## (Intercept)           x 
##       5.683       1.966 
## 
## Degrees of freedom: 100 total; 98 residual

abline(a = fit$coef[1], b = fit$coef[2], col = "red", lty = "dashed", 
    lwd = 3)

(fit <- rq(y ~ x, tau = 0.1))

## Call:
## rq(formula = y ~ x, tau = 0.1)
## 
## Coefficients:
## (Intercept)           x 
##       4.285       2.220 
## 
## Degrees of freedom: 100 total; 98 residual

abline(a = fit$coef[1], b = fit$coef[2], col = "blue", lty = "dashed", 
    lwd = 3)